Week 10 : Monte Carlo
The University of Sydney
Important
This presentation is based on the SOLES reveal.js Quarto template and is licensed under a Creative Commons Attribution 4.0 International License.
The question was what are the chances that a Canfield Solitaire laid out with 52 cards will come out successfully? After spending a lot of time trying to estimate them by pure combinatorial calculations, I wondered whether a more practical method than “abstract thinking” might not be to lay it out say one hundred times and simply observe and count the number of successful plays.
In this example, the mathematical statistics of the procedure is:
\begin{align*} \mathbb{E} [X] = \int_A t f(t) \, dt \approx \frac{1}{N} \sum_{i = 1}^N X_i \end{align*}
From last slide, we can write the expected value of a random variable X as a sum:
\begin{align*} \mathbb{E} [X] = \int_A t\cdot f(t)\, dt \approx \frac{1}{N} \sum_{i = 1}^N X_i \end{align*}
We can replace X with a function g(X). Then the previous equation becomes:
\begin{align*} \mathbb{E} [g(X)] = \int_A g(t)\cdot f(t)\, dt \approx \frac{1}{N} \sum_{i = 1}^N g(X_i) \end{align*}
Here, we assume to draw X_i \sim f
Goal is to evaluate the following definite integral
\begin{align*} \int_0^1 e^{-x^2/2}\, dx \approx \frac{1}{N} \sum_{i = 1}^N g(X_i). \end{align*}
Find a random variable such that the above holds. Need to be careful such that:
One example is:
R with runifHow do we simulate random variables that can have any probability distribution?
In R, we can draw random Gaussian variables using the function rnorm but how do these functions really work?
Two methods (others exist):
\begin{align*} F(x) = \int_{-\infty}^x f(t) \, dt \end{align*}
\begin{align*} X = F^{-1}(U) \end{align*}
The exponential distribution has a pdf of the form:
\begin{align*} f(x) = \begin{cases} \lambda e^{-\lambda x} & \quad \text{if } x > 0\\ 0 & \quad \text{if }x \le 0 \end{cases} \end{align*} and a cdf of the form: \begin{align*} F(x) = \begin{cases} 1 - e^{-\lambda x} & \quad \text{if } x > 0\\ 0 & \quad \text{if }x \le 0 \end{cases} \end{align*}
The inverse cdf is:
\begin{align*} F^{-1}(p) = -\frac{\log(1 - p)}{\lambda}, \quad 0<p<1 \end{align*}
To sample from the exponential distribution, we can do the following:
f(x) is a normal mixture distribution, difficult to sample directly
M = 2.5
g(x) \sim N(0, 1)
library(wavethresh) # for dclaw
# Simulate data
n <- 1000
x <- seq(from = -3, to = 3, length.out = n)
f <- dclaw(x) # target density
g <- dnorm(x) # proposal density (standard normal)
M <- 2.5 # envelope multiplier
Mg <- M * g
# Acceptance-Rejection Sampling
accepted <- numeric(0)
n_accept <- 0
while (n_accept < n) {
# Sample candidate from proposal distribution
x_candidate <- rnorm(1)
u <- runif(1)
# Compute acceptance probability
if (u <= dclaw(x_candidate) / (M * dnorm(x_candidate))) {
accepted <- c(accepted, x_candidate)
n_accept <- n_accept + 1
}
}
# Plot histogram of accepted samples
df_samples <- data.frame(x = accepted)
ggplot(df_samples, aes(x = x)) +
geom_histogram(aes(y = ..density..), bins = 40, fill = "skyblue", color = "white") +
stat_function(fun = dclaw, color = "red", size = 1.2) +
theme_minimal() +
ggtitle("Accepted Samples vs True Density (f(x))")Each Pop Mart Monster blind box costs $22.
There are 10 unique figures in the full set.
Each box contains one random figure.
❓ How many boxes would you need to purchase — and how much money would you expect to spend — to collect all 10 figures?
Expected number of shops you need to do:
\begin{align*} \mathbb{E} [T] = n \left(\frac{1}{1} + \frac{1}{2} + \cdots + \frac{1}{n}\right) \approx n \log n + 0.5772n + 0.5 \end{align*}
For the example, we have n=10 (the number of unique figures) so \mathbb{E}[T] \approx 30
On average, you need to buy about 30 boxes.
The expected total money spent is 30 \times 22 = 660 dollars
n.shops.sim <- c()
for(i in 1:1000) {
collected <- c(); n.shops <- 0
while(length(collected) < 10) {
newitem <- sample(10, 1)
collected <- union(collected, newitem)
n.shops <- n.shops + 1
}
n.shops.sim <- c(n.shops.sim, n.shops)
}
mean.sim <- mean(n.shops.sim)
mean.sim ## the average number of eligible shops.
ggplot(data.frame(x = n.shops.sim)) +
theme_minimal() +
geom_histogram(
aes(x = x, y = after_stat(density)),
colour = "black", fill = "white") +
labs(x = "Number of Blind Box", y = "Density") +
geom_vline(xintercept = mean.sim,
linetype = "dotted",
colour = "red")[1] 29.733
After 1000 draws, how many black balls do you expect to be in the urn?
onerun <- function(n = 1e3, return.last = FALSE) {
urn <- c("white", "black")
prop.blacks <- NULL
for(step in 1:n) {
flip <- sample(urn, 1)
if (flip == "white")
urn <- c(urn, "white" )
else
urn <- c(urn, "black" )
nblack <- sum(urn == "black")
nballs <- length(urn)
prop.blacks <- c(prop.blacks, nblack/nballs)
}
if (return.last)
prop.blacks[n]
else
prop.blacks
}
set.seed(5003)
m <- 5
n <- 1e3
df <- data.frame(x = 1:n,
polya.urn = unlist(lapply(rep(n, m), onerun)),
Simulation = rep(LETTERS[1:m], each = n))
ggplot(df) + theme_minimal() +
geom_line(aes(x = x, y = polya.urn, colour = Simulation)) +
labs(x = "Number of draws", "Proportion of black balls")An European call option gives you the right to buy a stock at a particular price at expiry
Example: 90 day call option of a stock with $105 strike. Let’s say a stock is currently priced at $100. If it hits $108 at day 90, then you can exercise the option and make a $3 profit from this stock. If its price is below $105 at day 90, then you make nothing.
Think of this as the insurance premium you’d pay to guarantee the right to buy later at a known price.
A put option gives you the right to sell a stock at a particular price at expiry
\text{Payoff} = \max(S_T - K, 0)
where S_T is the stock price at expiry and K is the strike price.
There are three possible scenarios:
Equation for the future price S of a stock:
\begin{align*} dS = S (\mu dT + \sigma \sqrt{dT} \mathcal N) \end{align*}
Simulate many possible future stock prices S_T based on the above equation.
For each simulation, calculate payoff.
The price you are willing to pay is the expected value of the payoff E(\text{Payoff}).
Example: 90 day call option of a stock with $105 strike. currently priced at $100.
set.seed(5003)
ndays <- 90
S0 <- 100 # Starting stock price
mu <- 0.05 # Stock drift rate per year
dT <- 1/365
sigma <- 0.5 # volatility
strike <- 105 # strike price of call option
# Draw some normal random numbers ~ N(0,1)
niters <- 5000 # Number of Monte Carlo simulations to do
option.payoff.arr <- c()
stock.prices <- list()
S90 <- c()
for(j in 1:niters) {
N = rnorm(ndays)
S = S0
S_arr <- c(S0)
for(i in 1:ndays) {
dS = S*(mu*dT + sigma*sqrt(dT)*N[i])
S = S + dS
S_arr <- c(S_arr, S)
}
S90 <- c(S90, S)
stock.prices[[j]] <- S_arr
# For a call option, if the stock price at expiry is
# higher than strike price, then we profit
# Otherwise, we don't exercise the option
option.payoff = max(S - strike,0)
option.payoff.arr <- c(option.payoff.arr, option.payoff)
}
mean(option.payoff.arr) ## expected price you are willing to pay for this option[1] 8.723601
# library(RColorBrewer)
m <- 5 # Number of stocks to plot
dat <- data.frame(Day = 0:ndays,
`Stock price (in $)` = unlist(stock.prices[1:5]),
Simulation = rep(LETTERS[1:m], each = ndays + 1),
check.names = FALSE)
ggplot(dat) + theme_minimal() +
geom_line(aes(x = Day, y = `Stock price (in $)`, colour = Simulation)) +
geom_vline(aes(xintercept = 90, clour = "blue"))